#=
测试势能的形状
=#


include("../src/pimary.jl")
using ..pimary
using Plots
using LinearAlgebra


function run1()
    xvals = collect(Base.OneTo(100))/25
    potes = [bare_Coulomb([x, 0.0, 0.0], [2.0, 0.0, 0.0]) for x in xvals]
    plt = plot(xvals/4, potes, xlims=(0.25, 0.55), ylims=(0.0, 20.0), label="Coulomb")
    pote2 = [Kelbg_repulse([x, 0.0, 0.0], [2.0, 0.0, 0.0]) for x in xvals]
    #println(pote2)
    plot!(plt, xvals/4, pote2, xlims=(0.25, 0.55), ylims=(0.0, 20.0), label="\$\\Delta\\beta = 0.1\$")
    #
    pote2 = [Kelbg_repulse([x, 0.0, 0.0], [2.0, 0.0, 0.0]; λ11=0.3162*pimary._WLElec) for x in xvals]
    plot!(plt, xvals/4, pote2, xlims=(0.25, 0.55), ylims=(0.0, 20.0), label="\$\\Delta\\beta = 0.01\$")
    savefig("test.png")
    #println(Kelbg_repulse([2.0, 0.0, 0.0], [2.0, 0.0, 0.0]))
end


function run2()
    Np = 3
    Nt = 1
    vol = 1
    fbup = WrdFb(Float64, Np, Nt)
    fbdn = WrdFb(Float64, Np, Nt)
    uppos = vol*rand(Np, 3)
    dnpos = vol*rand(Np, 3)
    for ip in Base.OneTo(Np); for it in Base.OneTo(Nt)
        fbup[ip][it, :] = uppos[ip, :]
        fbdn[ip][it, :] = dnpos[ip, :]
    end; end
    #
    for it in Base.OneTo(Nt)
        uppos = beads(fbup, it)
        dnpos = beads(fbdn, it)
        rpos = vcat(uppos, dnpos)
        phi = Kelbg_repulse_pairs(rpos)
        println("$(it) $(phi)")
    end
end


"""
测试对
"""
function run3()
    melec = 1.0
    Δβ = 0.10
    hbar = 1.0
    ther = sqrt(2π*(hbar^2)*Δβ/melec)
    t11 = sqrt((hbar^2)*Δβ/melec)
    #
    Np = 2
    Nt = 3
    vol = 1
    fbup = WrdFb(Float64, Np, Nt)
    uppos = vol*rand(Np, 3)
    #uppos[2, :] = uppos[1, :]
    for ip in Base.OneTo(Np)
        fbup[ip][1, :] = uppos[ip, :]
    end
    #每个it上运动一个随机的WL
    #offs = 0.5*ones(Nt, Np, 3)#rand(Nt, Np, 3)
    offs = rand(Nt, Np, 3)
    for it in 2:1:Nt
        for ip in Base.OneTo(Np)
            fbup[ip][it, :] = fbup[ip][1, :] .+ vol*(offs[it, ip, :] .- 0.5)*ther
        end
    end
    #每个it上的Kelbg势能大小
    println("ther: $(ther)")
    prevpot = zeros(Nt)
    for it in Base.OneTo(Nt)
        ri = fbup[1][it, :]
        rj = fbup[2][it, :]
        pot = Kelbg_repulse(ri, rj; λ11=t11)
        nij = norm(ri .- rj)
        println("$(it) $(ri) $(rj) $(nij)")
        println("$(it) $(pot)")
        prevpot[it] = Δβ*pot
    end
    #每个it上的Kelbg导数
    aderv = zeros(Nt)
    for it in Base.OneTo(Nt)
        ri = fbup[1][it, :]
        rj = fbup[2][it, :]
        r0i = fbup[1][1, :]
        r0j = fbup[2][1, :]
        pot = Kelbg_repulse_dbeta(ri.-rj, r0i.-r0j; λ11=t11)
        #throw(error())
        nij = norm(ri .- rj)
        println("$(ri) $(rj) $(r0i) $(r0j)")
        println("$(ri .- rj) $(r0i .- r0j)")
        a = ri .- rj
        b = r0i .- r0j
        println(a.-b)
        println("$(it) $(pot) $(ri.-rj .- r0i.-r0j)")
        aderv[it] = pot
    end
    #
    melec = 1.0
    Δβ = 0.101
    hbar = 1.0
    ther = sqrt(2π*(hbar^2)*Δβ/melec)
    t11 = sqrt((hbar^2)*Δβ/melec)
    #
    Np = 2
    Nt = 3
    vol = 1
    fbup = WrdFb(Float64, Np, Nt)
    #uppos = vol*rand(Np, 3)
    for ip in Base.OneTo(Np)
        fbup[ip][1, :] = uppos[ip, :]
    end
    #每个it上运动一个随机的WL
    #offs = rand(Nt, Np, 3)
    for it in 2:1:Nt
        for ip in Base.OneTo(Np)
            fbup[ip][it, :] = fbup[ip][1, :] .+ vol*(offs[it, ip, :] .- 0.5)*ther
        end
    end
    #每个it上的Kelbg势能大小
    println("ther: $(ther)")
    postpot = zeros(Nt)
    for it in Base.OneTo(Nt)
        ri = fbup[1][it, :]
        rj = fbup[2][it, :]
        pot = Kelbg_repulse(ri, rj; λ11=t11)
        nij = norm(ri .- rj)
        println("$(it) $(ri) $(rj) $(nij)")
        println("$(it) $(pot)")
        postpot[it] = Δβ*pot
    end
    println(postpot - prevpot)
    println(aderv)
end

run1()
run2()
run3()

